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Abstract. We study the positive recurrence of multi-dimensional birth-and- 
dcath processes describing the evolution of a large class of stochastic systems, 
a typical example being the randomly varying number of flow-level transfers 
in a telecommunication wire-line or wireless network. 

We first provide a generic method to construct a Lyapunov function when 
the drift can be extended to a smooth function on R , using an associated 
deterministic dynamical system. This approach gives an elementary proof of 
ergodicity without needing to establish the convergence of the scaled version 
of the process towards a fluid limit and then proving that the stability of the 
fluid limit implies the stability of the process. We also provide a counterpart 
result proving instability conditions. 

We then show how discontinuous drifts change the nature of the stabil- 
ity conditions and we provide generic sufficient stability conditions having a 
simple geometric interpretation. These conditions turn out to be necessary 
(outside a negligible set of the parameter space) for piece- wise constant drifts 
in dimension 2. 



1. Introduction 

We study the stochastic stability of multi-dimensional birth and death processes 
X = (Xi, . . . , Xn) on , (N being an integer greater than 1) with state-dependent 
birth and death rates (respectively X(x) = (Xi(x))i=i...N and (<pi(x))i—i,,,N, with 
x = (xi, . . . ,xn)) being 0-homogeneous functions, i.e. such that X(ax) — X(x) for 
any a > and for any x € 1+. 

This assumption arises naturally in many queueing networks representing com- 
munication or manufacturing systems where the server capacity depends on the 
states of all queues. Cellular radio networks are a typical example: available trans- 
mission rate for customers in a particular cell is decreasing when the number of 
customers in the neighboring cells increases [2]. More generally, many telecommu- 
nication networks can be represented (at a sufficiently large time scale) as processor 
sharing networks, with processing rates that may depend on the number of cus- 
tomers at each node of the network [3] . The notion of fairness between the classes 
of customers has led to the introduction of allocations of bandwidth naturally de- 
pending on the proportion of customers of each class rather than on the number of 
customers [15]. This justifies our assumption of 0- homogeneity. 

l 
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The assumptions made here are also relevant for describing load balancing schemes 
between a set of servers or computer systems. Very few results have been discussed 
in the literature on simple schemes such as joining the shortest queue when the 
death rates are not constant. (For the constant case, see [TO]). 

A general framework for analyzing stochastic stability consists in applying the 
Foster-Lyapunov criteria, which are based on finding a suitable test function having 
a positive or negative mean drift in almost all states of the state space [MJ [HI 
116] . When there are further restrictive assumptions made on service rates, an 
appropriate Lyapunov function has been found in many cases. This is, for example, 
possible for rates being the solution of specific optimization problems [3] , or for small 
dimensions (2 and 3) when the rates of the process are constant on sub-faces of the 
orthant 9 . For more complex systems, however, finding a Lyapunov function can 
be a formidable task. 

An alternative tool for deriving stability conditions is to study whether the sys- 
tem of interest is stochastically comparable to another system that is easier to 
analyze. This approach was first used in the multi-class queueing context by Rao 
and Ephremides [19] and Szpankowski |22j . and later refined by Szpankowski [23], 
to characterize the stability of buffered random access systems. It was later gen- 
eralized to birth-and-death processes with state-dependent transitions with fixed 
birth rates and decreasing death rates with uniform limits in |4J. This approach 
relies however on quite specific assumptions that are not verified even for simple 
processor sharing systems. 

Finally, many stability results have been obtained using the so-called ODE (ordi- 
nary differential equations) methods. A powerful exposition of these ideas applied 
to controlled random walks can be found in Chapter 10 of [18] and in [llj . The 
use of ODE methods is usually coupled with the analysis of fluid limits: first the 
convergence of a scaled version of the process towards a fluid limit is proven; then 
(under restrictive conditions) the stability of the fluid limit is proven to imply the 
stability (positive recurrence in our case) of the stochastic process. Stability condi- 
tions for a wide class of multi-class queueing networks with work-conserving service 
disciplines [6l [7] have been derived using these steps. In many recent papers (see, 
for instance, [32]), very involved proofs have been considered to demonstrate that 
the state of the stochastic network (under an appropriate space-time scaling) con- 
verges to a deterministic system whose evolution is represented by a differential 
equation of the form: 

±x(t)=S(x(t)), (1) 

where S is the drift of the stochastic process and x{t) corresponds to a fluid model. 
It is worth mentioning that the proof of such a convergence found in [12j holds for 
state processes that are not necessarily Markovian. 

It turns out however that such a convergence does not hold in general when 
the drift cannot be extended to a continuous function on Wf. When the drift 
vector-field is discontinuous, the trajectories of a fluid equivalent system enter slid- 
ing modes and the differential equation has to be replaced by a new dynamical 
system defined piece-wise by differential equations -^x(t) — S(x(t)), where 5 is a 
convex combination of drifts of points in the neighborhood of the discontinuity. 
Such a phenomenon was already emphasized in [9] where 5 was called the "second 
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vector-field" (see also PSE]). Unfortunately, 6 is difficult to compute in many 
cases (it depends crucially on the statistical assumptions made) and has not been 
characterized in general. 

Our contribution is three-fold. First, we present a way to construct a Lyapunov 
function using the deterministic differential equations driving the fluid limits dy- 
namics when the drift is continuous. This is done without proving that a determin- 
istic differential equation indeed represents the behavior of a scaled version of the 
underlying process. Such proofs have been necessary so far while our approach is 
more direct and elementary for proving stochastic stability and gives another sim- 
ple understanding of the meaning of such fluid limits for obtaining the stochastic 
stability. Moreover, it provides a systematical way of finding Lyapunov functions 
for stochastic systems. The advantage of finding a Lyapunov function explicitly is 
that it potentially gives much more precise information on the nature of the con- 
vergence of the process towards its stationary regime 16 . We also show that in 
the case of conservative drifts, the complexity of the problem can be considerably 
reduced. We then give a counterpart result, leading to instability. 

Second, we give general necessary conditions for stability in the case of discon- 
tinuous drifts. These conditions have a natural geometric interpretation. 

Third, we use these conditions to get a sharp geometric characterization of the 
stability set in case of piece- wise constant drifts in dimension 2. We give in partic- 
ular an algorithm allowing to conclude if the process is stable or not, for all fixed 
birth rate parameters outside a set of dimension 1. 

The paper is organized as follows. In Section [2] we describe the model in detail 
and discuss the methodology used in the subsequent analysis. In Section [3] we ex- 
amine the case when the drift vector-field can be extended to a continuous function. 
Section [4] is devoted to deriving sufficient stability conditions in the case of discon- 
tinuous drifts. We start by showing that fluid limits in this case are cumbersome 
and then proceed to presenting our approach in a generic scenario. In Section [5l 
we show how this approach may be applied to the processes in dimension 2 with 
piece-wise constant drifts in order to obtain a sharp geometric characterization of 
the stability region. Section [6] illustrates our various results and show that our 
sufficient conditions are not necessary in dimension 3. Section [7] concludes the 
paper. 

2. The model 

Let N be an integer greater than 1. We denote by the positive orthant of 
A (where A in this paper will be Z or K) while A+ * stands for Af \ {0}. 

Let ej be the vector of Z+ defined by (e,), = 1, = 0,3 ^ *• If n °t specified 

otherwise, | | denotes the usual Euclidian norm. The notation x < y is used for 
the coordinate- wise ordering: Vi, Xi < j/i, and we denote by (x,y) the usual scalar 
product of two vectors in M. N . A process X or a trajectory u started in x at time 
will be respectively denoted by X x and u x . 

Assume that AT is a continuous-time Markov process on Z^ with the following 
transition rates: 

q(x,x + e<) = Xi(x), 
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q(x,x ~ a) = (j>i(x), 

where A = (\i)i=i..,d and <ft = {4>i)i=i.~d are vectors consisting of positive 0- 
homogeneous functions from W\_ to R. 

The drift function 5 = (Si, ... , Sn) = \ — (f> is bounded, which guarantees that the 
process X is non-explosive. Hence we may assume that X and all other stochastic 
processes treated in the sequel have paths in the space D = D(R + ,Z^) of right- 
continuous functions from R + to with finite left limits. Recall that a stochastic 
process with paths in D can be viewed as a random element on the measurable space 
(D, T>), where V denotes the Borel c-algebra generated by the standard Skorokhod 
topology [H]. 

We are interested in conditions on the drift vector-field 8 ensuring that the 
process is either stable (recurrent) or unstable (transient or null- recurrent). In 
subsequent sections we shall find such conditions with the use of the so-called 
Foster-Lyapunov criterion (see, e.g. [24} l9l H6]). 

3. Smooth drift 

3.1. Stability conditions. Let (ODE) x be the following deterministic differential 
equation 

ju[t) = 8{u{t)). (2) 
u(Q) = x. 

We denote by u x a solution with an initial condition x. Dehne S = {x : \x\ = 1}. 

Theorem 1. Assume that S is a continuously differentiable function from R+ to 
Mr?. Assume that for allx G S, there exists a solution of (ODE) x such thatxit) = 
for all t > T x where T x < oo. Assume in addition that 

supTj; < oo. (3) 
xes 

Then X is positive recurrent and x i— > T x is a Lyapunov function. 

Proof. Note first that the homogeneity of the drift implies that if u x (t) is a solution 
of (ODE) x , then we can define u Kx (t) — Ku x (t/K) as a solution of (ODE)kx- 
Indeed: 

u Kx {0) = Kx, 

^-u Kx (t) = K^-u x (t/K) = S(u x (t/K)) = S(Ku x (t/K)) = S(u Kx (t)). 

Because 8 is C 1 on M.* , the flow (t, x) — > u x (t) is C 1 onRxRf , i.e. continuously 
differentiable in t and x (see for instance Theorem 1, page 299 in |13j). 

Define by F(x) = T x the time needed for u x to hit 0. We are going to show that 
F is a suitable Lyapunov function for proving positive recurrence of X. Due to the 
assumptions of the Theorem, F is a positive finite function. It is easy to see that 

^F(u x (t)) = -1 for all t < T x . 

Indeed, the difference between F(u x (t + h)) and F(u x (t)) is negative and is equal 
in absolute value to the time needed to reach u x (t + h) from u x (t), which is exactly 
h. Hence, the latter equality follows. 
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Moreover since u Kx (t) = Ku x (t/K), it follows that F{Kx) = KF(x), i.e. F 
is a 1-homogeneous function. This, in particular, implies that F(x) — > oo, when 
| a; | — > oo. The drift of F is given by 

AF(x) = J2Q(x,y)(F(y) - F(x)). 
y 

It will be convenient to approximate the variations of F using its gradient. We 
thus need to show that F is differentiable. For that purpose, define for each k > 
the functions 

T x>R =M{t: \u x (t)\ <k}. 

It is clear that T XjK — > T x for each fixed x as K — + 0. It is also clear that for each 
fixed x functions T XjK increase when k decreases. Note that, due to the continuity 
of u x (t) in t, 

\u x {T x , K )\ =K. 

Examine this equality: u x (t) is a differentiable function, while the norm of a differ- 
entiable function is also differentiable. Hence, we conclude that T x K is differentiable 
for all values of k > 0. In order to prove the differentiability of T x (and F(x)), it 
remains to show that the convergence T X<K — > T x is uniform in x. It follows from 
the following sequence of equalities: 

sup|T X)K -T x | = sup \T X \ — k- sup \T X \ — > 

X x:\x\=K x£S 

as k — > 0, due to ([3|). The last equality in the sequence follows from 1-homogeneity 
oiT x . 

As F is 1-homogeneous and differentiable, we have 

AF(x) < (VF(i), S(x)) + e, (4) 

for |a;| large enough, where e > may be made sufficiently small. 
Observing that 

M = (vF Wl )),^)HvWPWt))), 

we obtain, taking t — > 0, that 

(VFW.ifa)) < -e, 

For |x| large enough. Summarizing, the function F is such that F(x) — > oo as 
\x\ — > oo and that A_F(a;) < — e, for |x| large enough. After an examination of the 
Foster-Lyapunov criterion, it can be concluded that X is positive recurrent. 

□ 

Remark 3.1. Note that our approach is very similar to that of the fluid-limits 
approximation as it is developed for instance in |llj . However, we do not need to 
prove that the differential equation we look at is indeed the one that represents the 
behavior of our birth- and- death process on the fluid scale. Instead, we have clearly 
used that the differential equation (ODE) x is the characteristic equation of the PDF 
(partial differential equation) (\7F(x), 5(x)) — — e which directly leads us to finding 
a Lyapunov function. 
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Remark 3.2. It is worth mentioning that our analysis is valid thanks to the 0- 
homogeneity assumption. Without this assumption, the stability of the birth and 
death processes cannot be described using the differential equation |TS|). A counter- 
example ( with continuous drifts ) can be found in [4] . 

Remark 3.3. We believe that our construction could be considered for continuous 
drifts vector-fields using approximation arguments. This falls however out of the 
scope of this paper. 

3.2. Instability conditions. In this section, we consider a reverse statement es- 
tablishing instability relying on the previously considered dynamical system. It is 
much more challenging to state generic instability conditions based on the ODE (J2} , 
without a direct use of fluid limits i.e. without needing to prove the convergence of 
a scaled version of the process towards the trajectories of the ODE (see [T7] for the 
construction of a Lyapunov function proving the transience of multi-class queuing 
networks with routing). 

We use hereafter explicitly the convergence to the fluid limits and refer to [12] 
(see also in the discrete time setting [11]) for a proof of this result. The next theorem 
is hence essentially a combination of proving the convergence of the scaled process 
and the extended version of the (instability part of) Foster-Lyapunov criterion (see 
for instance Theorem 2.2.7 in [9]). 

Theorem 2. Assume that 5 is a Lipschitz function outside a neighborhood of the 
origin. Assume further that there exists a strictly positive time T and a number 
a > 1 such that for all x, with \x\ — 1, a solution u x of {ODE) x is defined on an 
interval [0,t x ] with t x > T and verifies \u x (T)\ > a. Finally suppose that for all 
times t and points x 



lim P(sup 

K — >oo T 



X Kx {Kt)) 



Kt 



> e) = 0. 



Then X is transient. 



Proof. Note first that the Lipschitz condition ensures that the ODE @ has a 
unique solution u x for each x £ » on an interval [0,t x ]. Furthermore, since is 
necessarily a stable point if it is an equilibrium point, the conditions of the theorem 
imply that the trajectories u x (t) did not hit before time T. Using [12 [Theorem 
4.1] (which actually needs only that 6 is continuous), we have the convergence of 
the process towards its fluid limit in the sense that 

X Kx {Kt) 



lim P 



K 



'■(t) 



>e =1, 



for each interval [0,t] included in the interval [0,t x ] where the ODE has a solution. 
This implies that for all x, there exists K x such that 

E\X K * X {K X T)\ - K x \x\ > (u x {T) - l)K x + e > 0. 

Hence 

E\X K * X {K X T)\ - K x \x\ >{a- l)K x . 

We now need to prove that sup K x = K < oo to be able to make use of the extended 
Foster-Lyapunov criterion (see for instance [9] [Theorem 2.2.7]). Denote by B v (x) 
a ball with radius r] and center x. 
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Fix r\ > 0. Using the martingale decomposition, X can be decomposed as 



x '" (t> = ¥4f i<xKl(s))d ' + 

i r Kt 

= x + — S(X Kx (s))ds + 
& Jo 



K 



M Kt 
K 



M 



Kt 



K ' 



where M is a martingale that satisfying: 

Tpf Mk * 

\0<s<t A 



< A 



K 



1/2 



< r/. 



for K large enough and A being a positive constant. Defining now = 
£[sup < s < T ±\X Kx {Ks) - X K y(Ks)\], we get that: 



f K {t) = \ x - y \ + 2rj+-E 







sup 


/ 5(X Kx (u))du 


0<s<T 


Jo 



< 3r] + E 


sup 


I 




.0<s<T 





X Kx {Ku) x K y(Ku) 



K 



-)-5(- 



K 



)du 



where we have successively used the convergence of the martingale part of the 
process to and the homogeneity of the drift. We can now condition on the 
event A e = {Vs;X Kx {s) <£ B € (0)} n {Vs : X R y(s) £ B e {0)} and use the Lipschitz 
condition. There exist two constants B and C such that: 



fx{t) < 3 V + BE [ 
Jo 



X Kx {Ku) X K y(Ku) 



K 



K 



du + CP{A t ). 



Using the theorem assumptions, we get that there exists K such that for all K > 
Kq, P(A e ) < V- We now conclude from Gronwall's lemma that there exists a 
constant D such that for all y G B v (x): 

M*) < D v, 

which in turn implies that supK x = K < oo. The last assertion concludes the 
proof. 

□ 



3.3. Gradient systems. Finding an explicit form for the function F might be dif- 
ficult in general. However, under slightly stronger assumptions on the drift function 
5, we can construct an explicit Lyapunov function directly from the vector-field 5. 
This is a well known fact in the theory of deterministic dynamical systems. 

Proposition 1. Assume that 5 is a conservative vector-field i.e., S — — W, and 
assume that 

V{x) > a > 0. (5) 
Assume further that there exists e > 0, such that 

\o~(x)\ > e for all x. (6) 
Then V{x) is a Lyapunov function and X is positive recurrent. 
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Proof. Note first that since 8 is O-homogeneous, V is 1-homogeneous. Using the 
1- homogeneity, there exists K < e 2 such that for |x| large enough, we can estimate 
the drift of V by: 

AV(x) < (VV(x),6(x)) +k 

= (V7(i),-VF(i)} + /! 
= -\V V(x)\ 2 + k < -e'. 

Furthermore, |V(x)| — * oo for |x| — * oo since it is a 1-homogeneous and strictly 
positive function. We can therefore apply the Foster-Lyapunov criterion. □ 

Remark 3.4. A vector- field 8 (xi,X2) = (<Ji(xi, xz), ^(xi, X2)) (on a completely 
connected set) is conservative if and only if -£^82{x) — ^7^1 ( x )- 



4. Discontinuous drifts 

4.1. Complexity of the fluid limits. So far we restricted ourselves to the case 
when the drift vector-field is continuous. The situation changes dramatically when 
this condition is dropped. When the drift vector-field is discontinuous, the trajec- 
tories of a fluid equivalent system near a point of discontinuity may enter " sliding 
modes" and the differential equation ([TJ) has to be replaced by a new dynamical 
system defined piece-wise by differential equations 4fx(t) — 8(x(t)), where 8 is a 
convex combination of drifts around neighborhoods of the discontinuities. 

Let us give a simple example of this phenomenon. Consider the following transi- 
tions with fixed birth rates Ai , A2 and death rates given by the following bandwidth 
allocation: 

4>i(x) = 1x2=0 + ail X2 >o, <h( x ) = lxi=o + a 2 l Xl >a- (7) 

Suppose Ai < a\ and A2 < (1 — pi) + aip\. This condition is known to be 
sufficient for stability of such a model in dimension 2, and it has been obtained 
through different methods (see, for instance, [3 [9], [4]). 

The following proposition proved in [20] characterizes the fluid limit of X. 



Proposition 2. The process — ^ — - converges in distribution when K — ^ 00 to- 
wards a process x(t) satisfying the differential equations: 

j t x{t) = 8(x(t)), for x(t) > 0, (8) 

j t xx(t) = 0, for xi(t)=0, (9) 

^-x 2 (t) = X 2 -(l- Pl )+a 2Pl , /orxi(i) = 0. (10) 



The stability condition is easily interpreted when considering the convergence 
to of the obtained fluid limit. This example shows however that even in a very 
simple case the fluid limit satisfies an equation different from ([1]). 

In the next subsection we develop an approach that allows us to find stability 
conditions in the case of discontinuous drift vector-fields without the use of fluid 
approximation. 
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4.2. Sufficient stability conditions. This section is devoted to identifying rather 
general conditions on the drift vector-field ensuring stability even in the presence of 
discontinuities for 5. These conditions lead to useful geometric stability conditions 
in dimension 2, which are discussed in Section [5] 

We start by considering a general vector-field of O-homogeneous drifts such that 
the number of discontinuities is finite. We construct a Lyapunov function by pasting 
together local Lyapunov functions and using a smoothing technique. This method 
was first used in [5]. 

We need to introduce a few more notations. Define by B £ (x) a closed sphere 
with radius e and center x. For a point x, denote by T> £ (x) the set of drifts in a 
neighborhood of x, i.e., 

V s (x) = {S(y) : y e B e (x)}. 
We then define by V* e (x) the set of vectors 

K,e( x ) = & G Rd : < ~ a ' W G T> £ (x)U{~x}} . 

We now state an assumption on the vector-field S(x) that we shall prove to be 
sufficient to characterize the stability region of the process. 

Assumption [Ax): For all x ^ 0, there exists e > and a > such that: 

KM + - 

The theorem below is the main result of this section. 

Theorem 3. Assumption (Ai) implies that X is positive recurrent. 

Before presenting a rigorous proof, we would like to explain the result intuitively. 
If the sets T> e {x) are finite for all x, assumption (Al) may be better understood 
using a simple geometric interpretation. Using Farkas' lemma, we can state that 
either assumption (Ai) is true or x is in the cone induced by the vectors of T> £ (x), 
i.e. there exist non-negative weights such that 

a. t S l = x. 

iei 

It is hence natural to expect that if x is never contained in the cone induced by 
the drifts S(y) at points y close to x (which is exactly assumption (Al)), then the 
process is stable. 
Proof. 

Consider first the function H and the vector-field v u constructed by 

Vu(x) = arg max (x,r]). 

H(x,u) = max (x,rj). 

The function H is a natural candidate for a Lyapunov function but it is discon- 
tinuous which complicates drastically the drift calculations and precludes having a 
negative drift in all points. We overcome this difficulty by considering a smoothed 
version of H. Let K e be a C°°-probability density supported on the sphere B £ (0) 
and introduce 

F{x) = / H(x,u)K e (u)du. 

JueB{e) 
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Then F is clearly a C°°-function. We shall now prove that F is a suitable Lyapunov 
function. First, notice that by assumption (Ai), H(x,u) > a, and hence, F(x) > a 
for all x. This, together with the observation that the function F is 1-homogeneous, 
implies that F(x) — > oo as \x\ — > oo. 

Note further that the compactness of the sphere guarantees the existence of u < e 
(arbitrarily small) such that for all x, 

\VF{x)-v u (x)\<-e. (11) 

Using assumption (A±) again, we get that 

(S{x),v u {x)) < -a. (12) 

Hence using fTTT 

(S(x),VF(x)) < (5(x), v u (x)) + Ce. (13) 

Again using the fact that F is 1-homogeneous, and combining [T2l and [T3l we get 
that 

AF(x) = V)(HV) ~ F (x)) < (<5(x), VF(x)} + e', 

y 

for large |a;|. Hence, F is a Lyapunov function. 

□ 

5. Piece-wise constant drift in dimension 2 

In this subsection we apply the general result of Theorem [3] to a particular case 
of a discontinuous drift function. We assume that the state space of the underlying 
process is N 2 and that the rate functions are piece- wise constant. Together with 
the assumption that the rate functions are O-homogeneous, this means that there 
is a finite number of cones where rate functions are constant. 

We start be introducing some notations that will be used throughout this section. 
Assume that there are N vectors vi, ...,Ujv such that v\ — ei,t>2 = &2 (where ei 
and e2 are vectors co-directed with one of the axes) and such that S(x) = S k for 
any x = Avu + Bvk+i- This means that the drift at any point of the cone defined 
by Wfe and Vk+i is equal to 6 h . 

Note that we do not require the vectors to be different. (This means that a 
cone reduces to a line when two consecutive vectors Vk and v k+i are equal.) 

Introduce also certain sets that will be crucial for the definition of the stability 
region. For each k — 1, N — 1, let 

Ul = {5 : 5 k = Av k + Bv k+l for some A > 0, B > 0, A + B > 0} 

and for each k = 2, N, let 

Ul = {5: aS k + (1 - a)5 k - 1 = Av k for some A > 0, a € [0, 1]} . 

We are ready to state the main result of this section. 

Theorem 4. Assume that 

S e SS = (u!f)u!) . 

k=l 

Then the Markov process X is positive recurrent. Conversely, if S belongs to the 
interior of the complement of SS , then X is transient or null-recurrent. 
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Proof. 

We start by proving the first part of the theorem. In order to do this, let us 
verify the conditions of Theorem[3J It is clear that only two situations are possible: 

(i) the vector x belongs to the interior of a cone defined by vectors Vk and Vk+i- 
In this case x — Avk + Bvk+i for some A > 0, B > 0. 

(ii) the vector x is collinear to a vector Vk- In this case x = Avk for some A > 0. 
Consider these two situations separately. In case (i), thanks to Theorem [3j we 

need to show the existence of a vector r\ such that 

{r),Av k +Bv k+ i) > 

and 

(r,,S k )<0. 

In geometric terms this means that there exists a vector rj such that vectors Avk + 
Bvk+i and S h belong to different half-planes separated by the line normal to vector 
r\. It is easy to see that the existence of such a vector is guaranteed by the fact that 

Consider now situation (ii). Applying again Theorem [3l we see that one needs 
to show the existence of a vector r\ such that 

{f],Av k ) > 0, 

and 

i^s'- 1 ) <0. 

If we interpret this again in geometric terms, it is equivalent to the existence of a 
vector r) such that vectors —Avk, S k and <5 fc_1 belong to the same half-plane defined 
by the line normal to r\. Direct computations show that it follows from the fact that 
5 U%. The proof of the positive recurrence under the assumption that S S SS is 
now complete. 

Let us now show that if 5 € int ^U/bJi 1 ^1 U ^l) ) then the Markov process X 
is not positive recurrent. We are going to prove that under the given conditions, 
the process is actually not rate stable which prevents stability. Assume that the 
process is started in a cone k and that 8 € U^. The strong law of large numbers 
(SLLN) then implies that with a positive probability the process stays in the cone 
{Av k +Bv k+1 ,A>0,B>0}. 

Assume that X is positive recurrent. This implies that ► 0, almost surely. 

Using the Martingale decomposition of the process, 

Using the boundedness of the transitions, the martingale Mt is such that E[Mf] < 
Ct, which implies the convergence in L 2 and in probability of ^ to 0, which in 
turns implies the almost sure convergence along a subsequence. Conditioning on the 
fact that the process stays in the cone we obtain j- f Q n 5(X x (s))ds — > 0, i„ — > oo, 
which combined with the ergodic theorem for the positive recurrent Markov process 
X implies that there exists a > such that 

= o^- 1 + (1 - a)5 k , 
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Figure 1 . Generic construction of the stability set for piece- wise 
constant drifts in dimension 2. 

which contradicts 5 £ U^. 

Suppose now that S £ U%. In this case the SLLN implies that with a strictly 
positive probability, the process stays in the set {Avk-i + Bv k ,A > 0,B > 0} U 
{Avk + Bvk+i,A > 0, B > 0}. Proceeding similarly as in the previous case, there 
exists a such that: 

= &5 k + (1 - &)5 k+1 , 

which contradicts 5 G U%. 

□ 

5.1. Fluid limits. A very interesting situation occurs within the framework of this 
section when there exist k and a £ (0, 1) such that aS k + (1 — a)5 k+1 — Av>~ for 
some A ^ 0. In this case we know that the fluid limit^3 with an initial state in 
the cone defined by Vk~i and or in the cone defined by and Vk+i enter a 
so-called "sliding mode" with their trajectory reaching the ray defined by Vk after 
a finite time and not leaving it after this time. The new drift 5 obtained during the 
sliding mode on Vk must be a convex combination of S k and S k+1 and must also 
be collinear with Vk- Hence, we can explicitly calculate S by solving the following 
system in a and A: 

aS k + (1 - a)<5 fe+1 = Av k . (14) 

The existence of a solution with a strictly negative A is necessary to get a stable 
system while the existence of a solution with strictly positive A is sufficient to get 
instability of the process, which corresponds respectively to the case where the fluid 
limits converge to zero or infinity. 

This is a very particular scenario, as in the case of a dimension higher than 
2, equation (|14|) is generally underdetermined and stability conditions cannot be 
characterized directly. 



for the existence of fluid limits, we refer to |20l and to |lll Theorem 1.2, in discret time 
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5.2. Algorithm. Recall that we defined the drift vector 6 to be equal to X(x) — 
4>{x). It is often the case in queueing and telecommunications applications that 
the function A (representing the arrival rate) is assumed to be a constant and the 
question is for which values of A the system under consideration is going to be stable. 
We present here an algorithm to construct the stability set when A is fixed. Define 
(within the description of the algorithm and the examples later on) 5 k — A — 4>k ■ 
The algorithm is given below. 

• step 1: Draw the points representing various values of ipk- 

• step 2 : Connect ipi to ip2, ip2 to ip^, etc. 

• step 3: For each k, draw the cone defined by vectors Vk and Vk+i based on 
point ip k - 

The compact set obtained is the stability region. 
A generic illustration of the previous algorithm is given in Figure [T] 



This section contains a number of examples illustrating the use of our results in 
various settings. 

6.1. Continuous drifts. 

Example 1. Consider a bandwidth sharing network representing a data commu- 
nication network as described in [3] with the following constraints on routing and 
capacity: 



where A and C are two convex sets. 

The following type of policies has been considered for many models in perfor- 
mance analysis (see for instance 1S\). Assume that the vector of traffic intensities 
and the bandwidth allocation are chosen such that 



where 5§ is the Minkowski function associated with a convex set S, i.e., 

Sg(x) — max(a;, it), 

and is also the Fenchel-Legendre transform of the support function l u es- It is now 
a matter of routine arguments to show that conditions [5| and [6| are satisfied if the 
interior of A(~\C is not the empty set. 

The results obtained in Section [3] allow to study numerically the positive re- 
currence of processes even when the drifts are too complicated to get an explicit 
solution for the associated ODE. 

Example 2. An example of a wireless network with two types of users competing 
for the same bandwidth could lead to the following deaths rates (using Shannon's 
formula and a state- dependent allocation policy): 



6. Examples 



A e A, 



x) — arg max (x, A — n) = — V<5*i 
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Figure 2. Example 3 : Drift vector fields for birth rates (0.4,0.8) 
and (0.5,0.8). 



where N is the thermal noise. 

Let us consider two possible vectors of arrival (birth) rate (0.4, 0.8) and (0.5, 0.8). 
The associated (ODE) can be solved numerically for any value Ai,A2 allowing to 
conclude for the positive recurrence of the process in the first case and the tran- 
sience in the second case, as shown in figure^ using the following properties of the 
trajectories of the ODE: 

• In the first case, all trajectories started from any state on the sphere hit 
in a bounded time, 

• In the second case, all trajectories started from the sphere do no reach 
a sufficiently small neighborhood of 0, from which we can conclude that 
all fluid-limits solutions stay outside of a ball of radius e and center 0. 
Moreover all trajectories do reach a state of norm bigger than 1 before a 
finite time T. 

6.2. Discontinuous drifts. We will give here a few examples where Theorem |4] 
provides interesting results. 

Example 3. We describe here how Theorem^can be used to obtain the well-known 
stability results for the so-called coupled-processors problem. Consider the allocation 
described as the most basic example with discontinuous drifts in Section^' 

4>i(x) = 1x2=0 + ail* 3 >ci, <h( x ) = lxi=o + a 2 l Xl>0 . 

It is clear that in this case the algorithm of Subsection \5.2\ allows us to recover 
the well-known stability region for this problem [5j [9] (see Figure E|): Ai < a\ and 
A 2 < (1 - pi) + a 2 pi or A 2 < a 2 and Ai < (1 - p 2 ) + aipi, with pi = 

Remark 6.1. Based on the previous example, one may think at first sight that 
the stability region is the smallest convex set containing all the vectors 4>{x) for x 
describing the state space. However, this is clearly not the case as it is illustrated 
in the generic example of Figure [7J In fact, the stability set may not be convex. 
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Figure 3. Stability region (set of birth rates A) for the 2 coupled 
processors of example 4 

Example 4. We now look at the same type of death rates as in the previous example 
but with different birth rates. Consider for instance a queuing or manufacturing 
system with 2 processors and two types of traffic: 

• some dedicated traffic arriving with intensity Xi to processor i, 

• some flexible traffic with intensity v that can be routed to either processor 
1 or 2 depending on the congestion level of both processors. 

Assume that the flexible traffic is actually routed to the processor with the smallest 
number of jobs in process (and to processor 1, say, if the processors are equally 
loaded, this last assumption having no impact on the stability conditions). Assume 
further that the arrivals of jobs of each type of traffic are following a Poisson process 
(independent of each other and of everything else) and that the processing times are 
exponentially distributed. The presence of interference or switching costs between 
different type of tasks raise the allocation of service (or death rates) described in 
the previous example with Oj < 1 . 

Using the notations of the previous section, the vectors Vk are: 

«i = (1, 0), v 2 = (1, 0), v 3 = (1, 1), Vi = (0, 1), v 5 = (0, 1), 

and the drifts are: 

S 1 = (\i—l,\2+v),6 2 = (Ai-ai, \ 2 +v-a 2 ),6 3 = (Xi+u—ai, A 2 -a 2 ), <5 4 = (Xi+u, A 2 

Using the results of the previous section, the interior of the stability region can 
be described as follows. 

Assume first that Xx — a\ > v + X 2 — a 2 - Then one of the two following conditions 
should hold: 

• X 2 + v < a 2 and Ai < 1 + ^7 i (A 2 + v), 

• or X 2 + v > a 2 + Ai — ai and A 2 — a 2 < v. 
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Symmetric conditions with the indices 1 and 2 reversed should hold if Ai — <zi + v > 

X 2 — CL 2 . 

6.3. Bounds on the stability region of 3 coupled processors. Consider a 
process of dimension 3 where the death rates of each dimension depend on whether 
the other coordinates are strictly positive or zero, so that for all i 7^ j ^ k, Xi > 0: 



4>i(x) = < 



a u Xj = 0, x k = 0, 
dij , Xj > 0, Xk 0, 
1, Xj > 0, Xk > 0, 



which leads to the following drifts: 

!S l : 8\ = Xi- ai, 5) = X 3 , 8\. = X k , for Xj = 0, x k = 0, 
S lj : Sf = Xi - dij, Sf = Xj - dji, 6% = Xk, for Xj > 0, x k = 0, 
5 = (Xi - l)j=i... 3 for x 3 > 0, x k > 0. 

Let us assume > > 1, so that <fi = (4>i, 4>2, ^3) is partially decreasing. 

Theorem 4.4.4 in [9] and Theorem 3 in [4] show that the stability region is a union 
of six regions corresponding to the six possible permutations of the coordinates. 
The first of these regions corresponding to the identity permutation is the set of 
(Ai, A2, A3) such that 

Si < 0, (15) 
Sf < Ai(l - a 23 ), (16) 
^ttoo + S1 3 tt 1q + 5f ttoi + <W < 0, (17) 

where 

^00 = P(Yi = 0,Y 2 = 0), TToi = P(Yi =Q,Y 2 > 0), 
7r 10 = P(Yi >0,Y 2 = 0), mi = P(Y 1 >0,Y 2 > 0), 

and Y = (Yi, Y 2 ) is a random vector distributed according to the stationary number 
of a process in which coordinate 3 would be always strictly positive. 

On the other hand the sufficient conditions obtained in Section [4] can be written 
as the complement of the following set: 



<5 > 0, (18) 

or 3i, j, and ai,a 2 > such that (aiS tJ + a 2 6, + e ) > 0, (19) 

or 3i and (aj)j=i...4 > such that 5jai + S'f a 2 + Sfa 3 + < 0. (20) 

It is not difficult to verify that conditions (jT&HTTj) are implied by the sufficient 
conditions (fT^l [2H)l obtained in theorem [3] while the converse is not true. 

7. Conclusions 

We derived various computable criteria of stability and instability for continuous 
drifts in any dimension and for piece- wise constant drifts in dimension 2, together 
with generic sufficient conditions for discontinuous drifts in any dimension. An 
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important direction of future research is to systematically characterize the second 
vector field and the stability conditions in dimension 3 and more. 
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